(loading matrix.ps...)print

% Coordinate System and Matrix Operators

% -  matrix  matrix
% create identity matrix
/matrix { 6 array identmatrix } def

% -  initmatrix  -
% set CTM to device default
/initmatrix {
    (initmatrix)=
    matrix defaultmatrix setmatrix
} def

% matrix  identmatrix  matrix
% fill matrix with identity transform
/identmatrix {
    dup type/arraytype ne { /identmatrix cvx /typecheck signalerror } if
    { 1 0 0 1 0 0 } exch copy
} def

% matrix  defaultmatrix  matrix
% fill matrix with device default matrix
/defaultmatrix {
    dup type/arraytype ne { /defaultmatrix cvx /typecheck signalerror } if
    DEVICE /defaultmatrix get exch copy
} def

% matrix  currentmatrix  matrix
% fill matrix with CTM
/currentmatrix {
    dup type/arraytype ne { /currentmatrix cvx /typecheck signalerror } if
    graphicsdict /currgstate get /currmatrix get exch copy
} def

% matrix  setmatrix  -
% replace CTM by matrix
/setmatrix {
    dup type/arraytype ne { /setmatrix cvx /typecheck signalerror } if
    graphicsdict /currgstate get /currmatrix get copy
    %pop
    (setmatrix)=
    ==
} def

%        tx ty  translate  -       % translate userspace by (tx, ty)
% tx ty matrix  translate  matrix  % define translation by (tx, ty)
/translate {
    dup type/arraytype ne
        { true 3 1 roll matrix } % no array: create array, concat later
        { false 4 1 roll } % array: do not create, do not concat later
    ifelse % bool tx ty matrix
    dup 0 1 put
    dup 1 0 put
    dup 2 0 put
    dup 3 1 put
    dup 4 5 4 roll put
    dup 5 4 3 roll put
    exch { concat } if
} def

%        sx sy  scale  -       % scale user space by (sx, sy)
% sx sy matrix  scale  matrix  % define scaling by (sx, sy)
/scale {
    dup type/arraytype ne
        { true 3 1 roll matrix }
        { false 4 1 roll }
    ifelse % bool sx sy matrix
    dup 0        % b sx sy mat mat 0
          5 4 roll put % b sy mat
    dup 1 0 put
    dup 2 0 put
    dup 3 4 3 roll put % b mat
    dup 4 0 put
    dup 5 0 put
    exch { concat } if
} def

%        angle  rotate  -       % rotate user space by angle
% angle matrix  rotate  matrix  % define rotation by angle degrees
/rotate {
    dup type /arraytype ne
        { true exch matrix }
        { false 3 1 roll }
    ifelse          % bool ang mat
    dup 0                 % b ang mat mat 0
          3 index cos put % b ang mat
    dup 1 3 index sin put
    dup 2 3 index sin neg put
    dup 3 4 3 roll cos put
    dup 4 0 put
    dup 5 0 put
    exch { concat } if
} def

% matrix  concat  -
% replace CTM by matrix \times CTM
/concat {
    graphicsdict /currgstate get /currmatrix get dup concatmatrix pop
} def

% matrix1 matrix2 matrix3  concatmatrix  matrix3
% fill matrix3 with matrix1 \times matrix2
/concatmatrix {
    3 copy type/arraytype ne exch type/arraytype ne or exch type/arraytype ne or
        { /concatmatrix cvx /typecheck signalerror } if
12 dict begin
    /D exch def
    aload pop 7 6 roll aload pop
    {l k j i h g f e d c b a}{exch def}forall
    a g mul c h mul add
    b g mul d h mul add
    a i mul c j mul add
    b i mul d j mul add
    a k mul c l mul add e add
    b k mul d l mul add f add
    D astore
end
} def

%        x y  transform  x' y'  % transform (x,y) by CTM
% x y matrix  transform  x' y'  % transform (x,y) by matrix
/transform {
12 dict begin
    dup type/arraytype ne { graphicsdict /currgstate get /currmatrix get } if
    aload pop % x y a b c d e f
    %(transform)=
    %pstack()=
    {f e d c b a y x}{exch def}forall
    a x mul c y mul add e add
    b x mul d y mul add f add
end
} def

%        dx dy  dtransform  dx' dy'  % transform distance (dx,dy) by CTM
% dx dy matrix  dtransform  dx' dy'  % transform distance (dx,dy) by matrix
/dtransform {
12 dict begin
    dup type/arraytype ne { graphicsdict /currgstate get /currmatrix get } if
    aload pop % x y a b c d e f
    %(dtransform)=
    %pstack()=
    pop pop % x y a b c d
    {d c b a y x}{exch def}forall
    a x mul c y mul add
    b x mul d y mul add
end
} def

%        x' y'  itransform  x y  % inverse transform (x',y') by CTM
% x' y' matrix  itransform  x y  % inverse transform (x',y') by matrix
/itransform {
12 dict begin
    dup type/arraytype ne { graphicsdict /currgstate get /currmatrix get } if
    %graphicsdict /currgstate get /scratchmatrix get copy invertmatrix transform
    aload pop {f e d c b a y x}{exch def}forall
    a d mul b c mul sub
    %(discr:)print dup =
    dup 0 eq { pop x y end /itransform cvx /undefinedresult signalerror } if
    1 exch div
    /invdet exch def
    d x mul b y mul sub c f mul add d e mul sub invdet mul
    c x mul neg a y mul add b e mul add a f mul sub invdet mul
end
} def

%        dx' dy'  idtransform  dx dy  % inverse transform distance (dx',dy') by CTM
% dx' dy' matrix  idtransform  dx dy  % inverse transform distance (dx',dy') by matrix
/idtransform {
12 dict begin
    dup type/arraytype ne { graphicsdict /currgstate get /currmatrix get } if
    aload pop % x y a b c d e f
    pop pop % x y a b c d
    {d c b a y x}{exch def}forall
    a d mul b c mul sub
    dup 0 eq { end /idtransform cvx /undefinedresult signalerror } if
    1 exch div
    /invdet exch def
    d x mul b y mul sub invdet mul
    c x mul neg a y mul add invdet mul
end
} def

% matrix1 matrix2  invertmatrix  matrix2
% fill matrix2 with inverse of matrix1
% revised using example at
%   https://groups.google.com/d/msg/comp.lang.postscript/JuclKJjfNRQ/hs5OV6Widu8J
/invertmatrix {
    2 copy type/arraytype ne exch type/arraytype ne or
    { /invertmatrix cvx /typecheck signalerror } if
12 dict begin
    exch aload pop {f e d c b a}{exch def}forall
    a d mul b c mul sub
    dup 0 eq { A end /invertmatrix cvx /undefinedresult signalerror } if
    1 exch div
    /invdet exch def
    d invdet mul
    b neg invdet mul
    c neg invdet mul
    a invdet mul
    c f mul d e mul sub invdet mul
    b e mul a f mul sub invdet mul
    7 -1 roll
    astore
end
} def

(eof matrix.ps\n)print
